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Abstract 

To process nonlinear and non-stationary signals, an extreme-point symmetric 
mode decomposition (ESMD) method is developed. It can be seen as a new alter- 
nate of the well-known Hilbert-Huang transform (HHT) method which is widely 
used nowadays. There are two parts for it. The first part is the decomposition 
approach which yields a series of intrinsic mode functions (IMFs) together with an 
optimal adaptive global mean (AGM) curve, the second part is the direct inter- 
polating (DI) approach which yields instantaneous amplitudes and frequencies for 
the IMFs together with a time- varying energy. Relative to the HHT method it has 
five characteristics as follows: (1) Different from constructing 2 outer envelopes, 
its sifting process is implemented by the aid of 1, 2 or 3 inner interpolating curves; 
(2) It does not decompose the signal to the last trend curve with at most one 
extreme point, it optimizes the residual component to be an optimal AGM curve 
which possesses a certain number of extreme points; (3) Its symmetry concept has 
wider extension than the envelop symmetry; (4) Its definition of IMF is more gen- 
eral; (5) The Hilbert-spectral-analysis approach is substituted by the data-based 
DI one. This new approach can not only yield clear distribution for instantaneous 
frequency and amplitude but also reflect the time-variation of the total energy in 
a distinct way. Besides the decomposition, the ESMD method also offers a good 
adaptive approach for finding the optimal global mean fitting curve to a given 
data. It is superior to the commonly used least-square method and the running- 
mean approach. We note here that the type of interpolation is not an essential 
factor, the approach developed here is also suitable for envelop-symmetric scheme 
adopted by the HHT method. 
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1. Introduction 

After two years exploration, now we put forward our "extreme-point 
symmetric mode decomposition (ESMD)" method for nonlinear and non- 
stationary signal processing. It can be seen as a new alternate of the well- 
known "Hilbert-Huang transform (HHT)" method developed by Huang et 
al. (1998) which is widely used nowadays. There are two parts for HHT 
method. The first and the key part is called "empirical mode decomposition 
(EMD)" which yields a series of intrinsic mode functions (IMFs), the second 
part is called "Hilbert spectral analysis (HSA)" which yields meaningful in- 
stantaneous amplitude and frequency for these IMFs. The scheme of EMD is 
to make a mode symmetric about its upper and lower envelopes interpolated 
by the local maxima and minima points, respectively. Physically speaking, 
it is a very natural phenomenon that a high-frequency small wave rides on 
a low-frequency big wave and the small one is almost symmetric about its 
crest and trough relative to the big one. Enlightened by this fact we choose 
another approach to do the decomposition. Our scheme is to make a mode 
symmetric about its own maxima and minima points. Differing from con- 
structing 2 outer envelopes, our sifting process is implemented by the aid 
of 1, 2, 3 or more inner curves interpolated by the midpoints of the line 
segments connecting the local maxima and minima points. We note here 
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that this approach had ever been tried in vain by Huang et al. (1998) for 
the 1-interpolating-curve case. As analyzed in the third section there are 
indeed some shortcomings for this case. In fact, 2 or 3 interpolating curves 
can lead successful decompositions. Besides the interpolating approach, the 
ESMD method has also some new developments in the stoppage criteria, the 
trend function, the concepts of symmetry and periodicity, the calculation of 
instantaneous amplitude and frequency together with the definition of IMF. 
By the way, as an alternate of HSA, a "direct interpolating (DI)" approach is 
also developed here for instantaneous amplitude and frequency. So the term 
"ESMD method" here is not merely confined to the decomposition, it also 
includes the postprocessing of IMFs. In the following we give some reviews 
on the related topics. 

1.1. About the Stoppage Criteria 

How to choose the stoppage criteria is a puzzling problem for the EMD 
users, after all, different sifting times may result in different decompositions 
[Huang et al. (2003)]. As stated in the review paper given by Huang and 
Wu (2008), how to optimize the sifting process is still an open problem. It 
is known to all that too low times of sifting may lead poor symmetry to the 
IMFs and render inaccuracy to the analysis of instantaneous frequency by 
Hilbert transform. Empirical speaking, to obtain a better symmetry it needs 
more times sifting. But it is not always true. Regardless of true or false, 
large number of sifting is not recommended [Huang et al. (2003), Wu and 
Huang (2009, 2010), Wang et al. (2010)]. Just as worried by Huang et al. 
(2003), too many times of sifting probably obliterate the intrinsic amplitude 
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variations and render the results physically less meaningful. In fact, this 
worry comes from the appearance of equal-amplitude IMFs as argued by 
the later researchers Wang et al. (2010). Through mathematical proof they 
have got the result that the upper and lower envelopes (in form of cubic 
splines) of a rigid symmetric IMF with sparsely populated extreme points 
have to degenerate to a pair of symmetric straight lines. In addition, the 
subsequent proof given by Wu and Huang (2010) indicates that the sparse 
condition is not necessary. Though this theoretical result is very attractive, 
our recent study [Wang and Li (2012)] shows that it is out of reach for the 
actual sifting process. Our results also indicate that the symmetry degree 
of IMFs changes in an intermittent manner, that is, as the sifting times is 
added, the sustaining-modulation state and sudden-turn state will appear in 
turn. So the symmetry of IMF will become better for the case of sustaining- 
modulation and worse for the case of sudden-turn. In addition, it follows 
from sifting tests that, after a certain sifting times, the average frequency of 
each IMF maintains unchanged or changes in an oscillating manner. Hence, 
the frequency ratio for neighboring IMFs can not decrease to 1 and the 
conjecture given by Wu and Huang (2010) is cracked. By the way, one can 
also reconsider the dyadic filter bank property of EMD [Flandrin, Rilling 
and Goncalves (2004), Flandrin and Goncalves (2004), Flandrin, Goncalves 
and Rilling (2005), Wu and Huang (2005)] and the frequency decomposition 
problem [Rilling and Flandrin (2008), Wu, Flandrin and Daubechies (2011)] 
for this limit case. 

According to the summary [Wang et al. (2010), Wang and Li (2012)], 
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there are four types of stoppage criteria: (1) the Cauchy type [Huang et 
al. (1998), Huang and Wu (2008)]; (2) the mean curve type [Rilling et al. 
(2003), Wang and Li (2012)]; (3)the S-number type [Huang et al. (2003)]; (4) 
the fixed-sifting-times type [Wu and Huang (2009, 2010)]. Among all these 
choices, if the mode symmetry is the sole attention in the sifting process then 
the mean curve type criteria are preferable, after all, the symmetry degree 
of IMFs changes in an intermittent manner along the sifting times. In this 
sense, the fixed-sifting-times type of criteria are not preferred if there is no 
prejudge on the symmetry. With the help of mean curve type criterion here 
we develop an ensemble "optimal-sifting-times" one to do the decomposition. 

1.2. About the Global Mean Curve 

For a given signal, its frequency analysis should be done on the os- 
cillating part. Hence, to clear away the global mean curve is the first and 
foremost problem. We note that the total mean (mathematical expectation 
in statistics) is just the simplest form of global mean curve. As illustrated 
by Huang and Shen (2005), the classical Fourier transform method is suit- 
able for the linear and stationary signal processing; the widely-used Wavelet 
transform method is suitable for the linear and non-stationary signal process- 
ing. Particularly, when the signal has a large evolutionary trend, the direct 
usage of Fourier and Wavelet transforms may lead distortion to frequency 
analysis. To get rid of the global mean curve, the least-square method and 
the running-mean approach are usually used. The least-square method may 
provide an optimal fitting curve for a given data in the sense of least vari- 
ance. But it is awkward in application, after all, it requires a priori function 
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form. The running-mean approach assigns the weighted mean value of sev- 
eral points to their center one, it may provide a smooth global mean curve 
to a given data. But this approach lacks of theoretical base and different 
choices of weight coefficients may result in different curves with uncertain 
boundaries. 

Since the EMD method adopts an adaptive scheme, to some extent, the 
global mean curve in form of trend function (with at most one extreme point) 
can be well extracted. But this kind of global mean curve may miss the evolu- 
tionary trend due to the bending limit. To make up this defect it only requires 
an superimposing management on several last lower-frequency modes. Cer- 
tainly, how to determine the number of modes is a question. Moghtaderi and 
his collages [Moghtaderi, Borgnat and Flandrin (2011), Moghtaderi, Flandrin 
and Borgnat (2011)] have discussed it by using the energy-ratio approach. 
By the way, since this approach is involved in the ratio of zero-crossing num- 
bers (identical to that of frequency) between the neighboring IMFs which 
is sensitive to the sifting times [Wang and Li (2012)], whether the obtained 
global mean curve is an optimal one or not is still a question. In fact, rather 
than decomposing the signal to the last trend function and superimposing 
the lower-frequency modes in return, we can directly stop the decomposition 
in a middle course. 

Differing from the EMD method, our ESMD method does not decompose 
the signal to the last trend function with at most one extreme point, it 
permits the residual component possesses a certain number of extreme points. 
One advantage of this processing is that: This kind of residual component 
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can reflect the evolutionary trend of the whole data much better and it can 
be understood as an optimal "adaptive global mean (AGM)" curve; Another 
advantage is that: We can optimize the sifting times by optimizing this AGM 
curve in the least-square sense. By the way, this optimization process itself 
is of great value. It offers a good adaptive approach for finding the optimal 
global mean fitting curve for a given data. This adaptive approach is superior 
to the commonly used least-square method and the running-mean approach. 
1.3. About the Concepts of Symmetry and Periodicity 

The type of symmetry is directly related to understanding the con- 
cepts of periodicity. Different from the envelop-symmetry used by EMD 
method, we adopt an extreme-point symmetry. This difference leads us to 
rethink the fundamental concepts of periodicity. 

For a constant function or a monotone function its periodicity and fre- 
quency arguments are not meaningful. Only when a quantity varies in a 
periodic oscillating manner, the frequency can be understood as an oscillat- 
ing change rate. The typical oscillating function is A cos out which accords 
with the ideal periodic variation of the substance in nature. Here A and oj are 
called the amplitude and frequency. Yet it is not always the case, such as a 
familiar damp vibration. Though its frequency (determined by the material 
property) maintains unchanged, its amplitude decreases as time goes on due 
to the resistance of air. This phenomenon is very universal. Certainly, the 
vibrating amplitude may also increase if it gains energy. The functions with 
fixed frequency and varying amplitude in the form A{t)cosut are defined 
as "weighted-periodic function" in our previous works [Wang and Li (2006, 
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2007), Wang and Zhang (2006)]. Mathematically speaking, the concept of 
periodicity can be also enlarged to the general form A(t) cos6>(t), where 9(t) 
is a continuous function. For convenience, we call it "generalized periodic 
function". In fact, the aim of EMD sifting is to extract a series of IMFs 
of this form. Now that the IMFs are generalized periodic functions, they 
can be directly extracted by mathematical approach. In this way, Hou and 
his cooperators [Hou, Yan and Wu (2009), Hou and Shi (2011, 2012)] have 
already made some effective explorations. 

For a periodic function, since its amplitude and frequency are all fixed 
constants, its symmetric characteristic is very clear. It can be understood 
as envelop symmetry or extreme-point symmetry. But from the viewpoint 
of material movement, as a matter of fact, the oscillation occurs around 
the equilibrium position. So the extreme-point symmetry actually reflects 
the local symmetry about itself (all the midpoints of the line segments be- 
tween the local maxima and minima points lie on the zero line). For a 
weighted periodic function, its frequency is fixed but its amplitude changes. 
At this time, its local symmetric characteristic is unclear, though it may ap- 
pear envelop-symmetric on the whole. To understand it from the viewpoint 
of extreme-point symmetry, the equilibrium should also shift its location. 
Hence, 

A(t) coscot = [A r (t) + A e (t)} coswt, (1) 

in which A r (t) and A e (t) should be understood as the real oscillating am- 
plitude and the equilibrium shifting amplitude. That is to say, during the 
oscillating process the corresponding equilibrium also shifts its location in the 
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same frequency (see Fig. 2). In addition, A r (t) and A e (t) are not indepen- 
dent, after all, A e (t) cos out reflects the trajectory variation of the midpoint 
for A r {t) cosut. As for the generalized periodic function, not only the ampli- 
tude but also the frequency changes. At this time, its amplitude can be also 
understood as above. But the frequency can not be understood as a simple 
instantaneous one, after all, the shifting of equilibrium location may distort 
the real oscillating frequency. For this case, the total oscillation can be seen 
as a synthesis of two components: 

A{t)cos6{t) = A r (t) cos 9 r {t) + A e {t) cos 9 e {t). (2) 

Particularly, in case A e (t) = the real amplitude A r (t) would degenerate to 
a constant and the function becomes an equal- amplitude form A r cos 9 r (t). 

This understanding is helpful for revealing the underneath nonlinear mech- 
anism of a complex system. The shifting phenomenon of equilibrium location 
is probably caused by the interaction of vibrations with different frequencies. 
Corresponding to the mode decomposition, this can be reflected by the in- 
teraction of IMFs. This topic is an attractive one to be discussed. 

According to the number of interpolating curves, we classify the ESMD 
into ESMDJ, ESMD JI, ESMDJII, • • •. In fact, their difference lie in the 
request on the equilibrium variation A e (t) cos 9 e (t) which should be under- 
stood as an interpolating function with all the midpoints. ESMDJ adopts 
a rigid extreme-point symmetry in the sifting process which requires all the 
midpoints almost lying on the zero line, that is, A e (t) ~ 0. For this case, all 
the IMFs should almost degenerate to the equal-amplitude form A r cos 9 r {t). 
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From the viewpoint of physics, this strategy is too rigorous. ESMDJI extends 
the concept of extreme-point symmetry, which permits the location shifting 
of the midpoint in such a manner: The trajectory variation A e (t) cos 9 e (t) 
should be envelop- symmetric about its odd and even interpolating curves. By 
the way, these two envelops differ from the common positive and negative 
outer envelops, after all, they may change their signs alternately. The sifting 
tests show that this odd-even type of extreme-point symmetry for ESMDJI 
is almost equivalent to the outer envelop symmetry for EMD (see Fig. 17). 
ESMDJII gives further extension to the concept of extreme-point symme- 
try. It liberalizes the restriction on A e (t) cos 9 e (t) and only requires the sum 
of two interpolating curves be symmetric with the third one. Certainly, the 
restriction on A e (t) cos 9 e (t) can be also liberalized in this way with more 
interpolating curves. 

1.4- About the Instantaneous Frequency 

The definition of instantaneous frequency is a controversial issue [Huang 
et al. (2009a)]. Just as analyzed above, only when a quantity varies in a pe- 
riodic oscillating manner, the frequency can be understood as an oscillating 
change rate during the process of moving back and forth. So there is no 
local meaning for frequency at a given point. But as argued by Huang et 
al. (2009a) and the references therein, there is indeed frequency modulating 
phenomenon. Therefore, to accord with the generalized periodic function 
A(t) cos 9(t) in mathematics, the derivative form u(t) = dQ/dt is recom- 
mended. To be meaningful, it requires d9/dt > 0. However, for a decom- 
posed IMF with denotation x(t) this calculation is not trivial. To solve this 
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problem, Huang et al. suggested the Hilbert transform which is popular used 
nowadays: 

1 f°° x(t) 

y(t) = H[x(t)\ = -P / -ii-dt, (3) 

in which P indicates the principal value of the singular integral. With the 
Hilbert transform, the analytic signal is defined as 

z{t)=x{t)+iy{t)=A{t)e»«\ (4) 

where 

A(t) = y/x\t)+y\t), 9(t) = arctan (j^j . (5) 

Here A(t) is the instantaneous amplitude and 9(t) is the phase function. 
Furthermore, the instantaneous frequency can be calculated by u(t) = d6/dt. 

Essentially, the Eqn.(3) defines Hilbert transform as the convolution of 
x(t) with 1/t, therefore, it emphasizes the local properties of x(t). In Eqn.(4) 
the polar coordinate expression further clarifies the local nature of this repre- 
sentation: it is the best local fit of an amplitude and phase varying trigono- 
metric function to z(t) [Huang and Shen (2005)]. So, in this sense, Hilbert 
transform is superior to the Fourier transform, Wavelet transform and other 
analytical forms. However, this approach has a disadvantage in unsatisfying 
the quadrature request delimited by the well-known Bedrosian and Nuttall 
theorems. That is to say, to use Eqn.(4) y(t) should be a quadrature func- 
tion of x(t). In addition, to get a meaningful instantaneous frequency it also 
needs a hypothesis that the derivative of 6(t) exists. 

In fact, no matter how the integral transform is defined, it is actually a 
uniform running-mean processing. Now that the processing is done on the 
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data, why not calculate the instantaneous frequency from the data in a di- 
rect way? Historically speaking, there are only crude estimation methods 
which can not reflect the instantaneous change. Just as reviewed by Huang 
et al. (2009a), there is a fundamental "zero-crossing method" which has long 
been used to compute the mean period or frequency for narrow band signals. 
Of course, this approach is only meaningful for mono-component functions, 
where the numbers of zero- crossings and extreme points must be equal in the 
data. Huang et al. had generalized the zero-crossing method by improving 
the temporal resolution to a quarter wave period with a running-mean ap- 
proach. This improvement yields a better estimation for the frequency, but 
it is still incapable of reflecting the instantaneous change. 




-1.5 1 1 1 1 1 1 1 1 

1 2 3 4 5 6 7 

Time 

Figure 1: An example of periodic oscillation with intermittences. 

As an instantaneous frequency it should be capable of reflecting the inter- 
mittent case as in Fig.l, rather than excluding the adjacent-equal situation. 
While the oscillation maintains unchanged at the extrema, the zero-crossing 
or any other locations its instantaneous frequency should be all 0. Certainly, 
the frequency may lose its meaning on the junctures (belong to a null set). 
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It doesn't matter, after all, the data itself is not smooth. Now that the period 
should be defined relative to a segment of time and the frequency needs to 
be understood point by point, we can conciliate this conflict by interpolating 
method. With this understanding, we have developed a "direct interpolat- 
ing" approach for the calculation of instantaneous frequency which will be 
illustrated in a latter section. 

1.5. About the Definition of IMF 

The EMD method defines an intrinsic mode function (IMF) with the 
following two conditions [Huang et al. (1998), Huang and Shen (2005)]: 

(1) In the whole data set, the number of extreme points and the number of 
zero-crossings must either equal or differ at most by one. 

(2) At any point, the mean value of the envelope defined by the local maxima 
and the envelope defined by the local minima is zero. 

The first condition can be understood as: the IMF's local maxima and 
minima points are septal with no adjacent zero-crossings, all its maxima 
should be positive and all its minima should be negative. Just as stated by 
Huang et al. (1998), this request on oscillating manner is for the rationality 
of defining an instantaneous frequency. According to this restriction, the 
function given in Fig.l is not an IMF. But from the viewpoint of physics, 
the occurrence of this intermittence phenomenon is possible. Therefore, it is 
reasonable for liberalizing the restriction and counting it as an IMF. 

The second condition requires that an IMF must have envelope symmetry. 
This restriction is for convenience of deducing a meaningful instantaneous 
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frequency by means of Hilbert transform. Now that the Hilbert transform is 
abandoned and the "direct interpolating" approach is adopted, the restriction 
should be liberalized. In fact, just as stated in our previous paper [Wang and 
Li (2012)] the rigid envelope-symmetric IMFs are out of reach for an actual 
data processing test. Therefore, only if a decomposed component satisfies 
a certain error requirement it can be seen as an IMF. In addition, notice 
that the odd-even type of extreme-point symmetry for ESMDJI is almost 
equivalent to the envelop symmetry for EMD and the three-curve type for 
ESMDJII is more general than the envelop symmetry, the request can be 
also liberalized. 

Based on the above analysis, the definition of IMF can be extended with 
the following two conditions: 

(1) To count all the adjacent equal extreme points as one, the local maxima 
and minima points should be septal, all its maxima should be positive and all 
its minima should be negative. 

(2) It should be almost envelop- symmetric or extreme-point symmetric in the 
generalized sense. 

We note that the envelop symmetry and the odd-even type of extreme-point 
symmetry are very good types which yield IMFs with suitable amplitude 
and frequency modulations, and too low request on symmetry would lead 
difficulty to frequency and energy analysis. 
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2. Decomposition Algorithm for ESMD Method 

We begin to introduce the ESMD method with the decomposition al- 
gorithm. By the way, besides this one there is also an interpolation algorithm 
for frequency. Here only one-dimensional signals are considered. Certainly, 
there is a precondition for the processing that the sampling rate of the ob- 
servational instruments should be known. It is a common sense that the 
local maxima and minima points are septal with counting all the adjacent 
equal extreme points as one. For convenience of programming, in case there 
are several adjacent equal extreme points, we only choose the first one as a 
representative. Our program code is exploited on the Scilab platform and 
the algorithm is as follows: 

Step 1: Find all the local extreme points (maxima points plus minima 
points) of the data Y and numerate them by E; t with 1 < i < n. 

Step 2: Connect all the adjacent Ei with line segments and mark their 
midpoints by Fi with 1 < i < n — 1. 

Step 3: Add a left and a right boundary midpoints F and F n through a 
certain approach. 

Step 4: Construct p interpolating curves L-y, ■ ■ ■ , L p (p > 1) with all these 
n + 1 midpoints and calculate their mean value by L* = (Li + • • • + L p )/p. 

Step 5: Repeat the above four steps on Y—L* until \L*\ < e (e is a permitted 
error) or the sifting times attain a preset maximum number K. At this time, 
we get the first mode M\. 

Step 6: Repeat the above five steps on the residual Y — M\ and get M 2 , M 3 • • ■ 
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until the last residual R with no more than a certain number of extreme 
points. 

Step 7: Change the maximum number ifona finite integer interval [K min , K max ] 
and repeat the above six steps. Then calculate the variance a 2 of Y — R and 
plot a figure with a/aQ and K, here <To is the standard deviation of Y . 
Step 8: Find the number K which accords with minimum a/a on [K min , K max ]. 
Then use this K to repeat the previous six steps and output the whole modes. 
At this time, the last residual R is actually an optimal AGM curve. 

There are several questions associated with this algorithm. In the follow- 
ing we explain them one by one. 

According to the fourth step, we classify the ESMD into ESMDJ, ESMDJI, 
ESMDJII, • • ■. ESMDJ does the sifting process by using only 1 curve inter- 
polated by all the midpoints; ESMDJI does the sifting process by using 2 
curves interpolated by the odd and even midpoints, respectively; ESMDJII 
does the sifting process by using 3 curves interpolated by the midpoints nu- 
merated by 3k + 1, 3k + 2 and 3(k + 1) (k — 0, 1, • • •), respectively. Certainly, 
we can also define other schemes with more interpolating curves according 
to this method. 

In the fifth step, besides the permitted error e, we can also adjust the 
maximum sifting times K. On the one hand, if e is the unique controlling 
parameter, it may leads to endless loop to the decomposition; On the other 
hand, if K is the unique controlling parameter, we know nothing about the 
symmetric property of each mode. Perhaps a small number of sifting may 
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lead good symmetric to a mode. So the wise choice is to use them all. To 
obtain a series of relatively reliable modes, we can fix e to be a very small 
value and control the decomposing process by changing K on a finite integer 
interval such that the last residual R (AGM curve) is an optimal one. So 
Step 7 and 8 are very necessary. In fact, only when the fitting curve of the 
data is an optimal one, the remaind signal can be seen as an actual oscillation 
caused by a series of wave fluctuations. 

Denote the original data and the AGM curve by Y = {yi}f =1 and R = 
{ri\f =1 i respectively. Commonly, relative to its total mean Y = YliLiVi/N 
the variance of the data is defined as 



In the applications, we usually choose e = O.OOlcro and use the ratio v = g/gq 
to reflect the degree of optimization for the AGM relative to the common 
total mean. 

In addition, the third step is associated with boundary processing which 
is a "benevolent see benevolence" problem. In our program codes we have 
developed the linear interpolation method given by Wu and Huang (2009) 
and revised the interpolation styles for the too steep boundary case [see 
Appendix A]. This revision can make the boundary much more stable, even 
for the tests with 100, 000 sifting times [Wang and Li (2012)]. In the following 




(6) 



Here we define the variance relative to the AGM by 




(7) 
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we test the decomposition effectiveness according to ESMDJ, ESMDJI and 
ESMDJII and our attention is mainly focused on the second one. 

3. Performance of ESMDJ 

ESMDJ does the sifting process by using only 1 curve interpolated by 
all the midpoints of the line segments between the local maxima and minima 
points. Though this case has been discussed by Huang et a/.(1998), it is 
worth reemphasized under the viewpoint of ESMD. We test it by a simple 
example as follows. 

Example 1: Y(t) = e~°- u sin (nt/2 + tt/3) with < t < 20. 

This is a weighted-periodic function with fixed frequency and varying 
amplitude. A signal may maintain its oscillating frequency and increase (or 
decrease) its amplitude with gaining (or losing) energy. So it is a very natural 
thing to meet the weighted-periodic function in signal processing. A good 
sifting scheme is anticipated yielding a unique IMF with small decomposition 
error. 

The detailed constructing process of the interpolating curve is shown 
in Fig. 2. In case there is no variance-ratio investigation, we know nothing 
about the effect of sifting times. To try the decomposition with 20 times 
sifting, it yields a result in Fig. 3, which includes 4 modes and a residual R. 
There is a common feature for these modes that all their amplitudes almost 
maintain unchanging. In fact, it is due to the scheme of rigid extreme-point 
symmetry adopted by ESMDJ. Through a simple geometric proof we see 
a curve with equal amplitude in the form Asin9(t) (9(t) is an increasing 
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Figure 2: All the midpoints (circle 
points) of the line segments between the 
local maxima and minima points and 
their unique interpolating curve (bulk in- 
ner curve). 




Figure 3: The decomposition result for 
the weighted-periodic function given by 
ESMDJ with 20 sifting times, here the 
horizontal axis stands for the time (sec- 
ond). 



function) must be an extreme-point symmetric one. Vice versa, it is also 
true. Particularly, the Model is not only extreme-point symmetric but also 
periodic. In fact, it is an approximation of the function 0.6 sin (irt/2 + 7r/3) 
which carries the most periodicity of the original signal. This try shows that 
20 times sifting accords with the case with enough extreme-point symmetry 
which leads a slow efficiency to the decomposition. Now that 20 times sifting 
is too excessive, a lower times with lose symmetry request may be better, 
but there is no effective criterion for it under the EMD circumstance. One 
progress of our ESMD is on the adoption of variance ratio. 

Firstly, we plot the distribution figure of the variance ratio v = a/a 
along the sifting times (see Fig. 4) and find out the optimal sifting time 3 
which accords with the minimum value v = 99.8%; Secondly, we output 
the corresponding decomposition result with 3 sifting times. It follows from 
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Fig.5 that the result is better than that in Fig. 3 with 20 times sifting, where 
Model can be seen as an approximation of the original signal and the others 
can be seen as decomposition error. Certainly, this decomposition is still not 
perfect, after all, the error amplitude achieves 0.3. 



100,6 
100.fi 
I uu 4 
100 3- 
100.2 
100.1 
100.0- 
99.9 
99.8 




Sifting Time 

Figure 4: The distribution of variance 
ratio v = <j/<7q along the sifting times 
for the weighted-periodic function. 




Figure 5: The decomposition result for 
the weighted-periodic function given by 
ESMDJ with 3 sifting times, here the 
horizontal axis stands for the time (sec- 
ond). 



Example 2: A segment of wind data observed at sea with 20Hz sampling 
rate. 

It follows from Fig. 6 that 18 is the optimal sifting times in the interval 
[1,30]. The corresponding decomposition in Fig. 7 yields 12 IMFs together 
with a residual R with 40% variance ratio. It means the AGM is the best 
fitting curve of wind data on the interval [1, 30]. At this time, the IMFs still 
have amplitude-modulation phenomenon. If the sifting times is added, more 
and more equal-amplitude IMFs may appear. In consideration of physical 
meaning stated by Huang et al. (2003), we do not expect too many equal- 
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amplitude modes. So the lower sifting times with imperfect decomposition 
is preferred for ESMDJ. 
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Figure 6: The distribution of variance 
ratio v = u/cfq along the sifting times 
for the wind data. 





Figure 7: The decomposition result for 
the wind data given by ESMDJ with 
18 sifting times, here the horizontal axis 
stands for the time (second). 



In all, due to the adoption of a scheme with rigid extreme-point symmetry 
the decomposition efficiency of ESMDJ is not high. Besides the rigorous 
request of symmetry there is also another aspect. It follows from Fig. 2 that 
when all of the midpoints are used for interpolating the generated curve 
may contain almost the same number of extreme points as the original data 
which may subsequently enter into the second mode. This also leads low 
efficiency to the decomposition. By the way, though ESMDJ has this defect 
its AGM curve may be very good. Comparatively, the EMD method has 
a relatively high efficiency of decomposing. One reason is that, the request 
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of envelope symmetry is lower than that of rigid extreme-point symmetry; 
Another reason is that, the upper and lower envelopes are interpolated by 
almost half number of the signal's extreme points and their mean curve 
discounts the number by almost a half. In view of these facts, it is very 
natural to do the decomposition by using 2 interpolating curves. 

4. Performance of ESMDJI 

ESMD_II does the sifting process by using 2 curves interpolated by 
the odd and even midpoints, respectively. The detailed constructing process 
of the curves is shown in Fig. 8. At this time, the decomposition of Example 1 
is trivial since the original curve itself is a permitted mode. In the following 
we test ESMDJI by three examples and analyze its characteristics from three 
aspects. 
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Figure 8: All the midpoints (big circle points) of the line segments between the local 
maxima and minima points of the data (thin solid curve) and the mean curve (thick solid 
curve) of their odd and even interpolating curves. 

4-1. Decomposition Tests 

Example 3: Y{t) = - sin(87rt) + 1.5e" a2 * sin (1.9vrt + tt/20) + (t - 2) 2 , < 
t < 4. 
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This signal is composed of one periodic function, one weighted-periodic 
function and one parabola. In the following we do the decomposition with 
ESMDJI. From Fig. 9 and 10 we see the decomposition is perfect. Model 
accords with the periodic function; Mode2 accords with the weighted-periodic 
function; The residual R accords with the parabola. 
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Figure 9: The distribution of variance 
ratio v = u/cfq along the sifting times 
for the composed data. 
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Figure 10: The decomposition result for 
the composed data given by ESMDJI 
with 29 sifting times, here the horizontal 
axis stands for the time (second). 



In the following we re-decompose Example 2 with ESMDJI. From Fig. 11 
we see the variance ratio v = o jo§ attains its minimum value at 30, which 
can be seen as an optimal sifting times in the whole interval [1, 100]. Fig. 12 
shows the corresponding decomposition which is more distinct than that in 
Fig. 7 given by ESMDJ. The decomposed IMFs accord with the components 
of wind turbulence with average periods 3.8s, 1.5s, 0.6s,- ■ -. In addition, the 
last residual R is an optimal AGM curve, which reflects the fundamental 
evolutionary trend of the wind speed very well (see Fig. 13). Certainly, there 
is a necessity for us to compare with EMD method. So the test is also done on 
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the code eemd.m (downloaded from http:/ /rcada.ncu.edu.tw/class2009.htm) 



for the non-noise case. By the way, to be more objective, the default 10 times 
sifting given by Wu and Huang (2009) is substituted by 30 here. From Fig. 12 
and Fig. 14 we see the difference is very clear. This difference probably lies in 
the sifting scheme, the boundary processing and the programming. Certainly, 
it is difficult for us to judge which IMF set is more reasonable. But it follows 
from the residual comparison in Fig. 13 we see the global mean curve given 
by ESMDJI is better than the monotone form given by EMD method, after 
all, our one is found by optimizing approach. It also indicates that our 
AGM curve is almost equivalent to the sum of Mode5, Mode6, Mode7 and 
the last trend function of EMD. Generally speaking, the global mean curve 
is a component with maximum magnitude, its deviation would lead large 
distortion to the oscillating IMFs. So in this sense, ESMDJI is preferable. 
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Figure 11: The distribution of the vari- 
ance ratio v = er/ao along the sifting 
times for the wind data. 



Figure 12: The decomposition result for 
the wind data given by ESMDJI with 
30 sifting times, here the horizontal axis 
stands for the time (second). 
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Figure 13: Comparison between the 
AGM curve of ESMDJI (the curve with 
small circles), the trend function (thick 
solid curve) as well as the sum of R 
and Mode5-7 of EMD (thin solid curve), 
where the dotted curve stands for the 
wind data. 



Figure 14: The decomposition result for 
the wind data given by EMD method 
with 30 sifting times. 



25 



Example 4: The day-averaged air temperature data from May 10, 2008 to 
Nov. 3, 2011 downloaded from the National Climatic Data Center. 

From Fig. 15 we see the variance ratio has some septal stable intervals, 
such as [20, 29], [36, 40], [43, 48] and [72, 76]. On all these stable intervals the 
variance ratio is almost identical, which leads almost the same result to the 
decomposition. At this time the optimal sifting time is 29, even when the 
whole interval is prolonged to [1,200]. The decomposition in Fig. 16 shows a 
very good yearly evolutionary trend for the air temperature. In addition, 
the decomposed IMFs accord with the components with average periods 
66day, 35day, 17day,- • -, which can be understood as bimonthly, monthly, 
semimonthly - • ■ temperature oscillations. Further analysis can be done in 
the way given by Huang et al. (2009b) and Bao et al. (2011). 
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Figure 15: The distribution of variance 
ratio v — a/ao along the sifting times for 
the temperature data. 



Figure 16: The decomposition result for 
the temperature data given by ESMDJI 
with 29 sifting times, here the horizontal 
axis stands for the time (day). 
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4-2. Symmetry Characteristic of the Modes 

When the decomposed mode is a periodic function, such as Model in 
Fig. 10, its symmetry characteristic is very obvious. At this time, not only the 
mean value of the odd and even interpolating curves but also themselves all 
approximate the zero line. This is a true extreme-point symmetry. But when 
the decomposed mode is a weighted-periodic function, such as the Mode2 
in Fig. 10, its symmetry is an extreme-point symmetry of odd-even type, 
which only requires the symmetry between the odd and even interpolating 
curves. In fact, this phenomenon is very universal. For an actual signal 
component, not only its amplitude but also its frequency changes along the 
time, such as the Model in Fig.12. It follows from Fig. 17 that this case is 
almost equivalent to the envelope symmetry. Since the magnitudes of the 
odd and even interpolating curves are smaller than that of the upper and 
lower envelopes, the convergence speed of their mean curve to the zero line 
may be quicker. So the ESMDJI method may need less sifting times than 
that for the EMD method to reach a relatively stable state. 

4-3. Effect of Sifting Times to the Decomposition 

According to the ESMD algorithm Step 4 and 5, to add the sifting 
times may make the modes more and more symmetric until \L*\ < e. More- 
over, since the permitted error e is preestablished (such as e = 0.001er ), 
more times of sifting may imply more symmetric modes. So we can antici- 
pate a finite times of sifting such that all the modes satisfy this permitted 
error. This case accords with a relatively stable variance ratio v = o~/o~q. 
For Fig. 9 the stable interval is approximately [25,34]; For Fig. 11 and Fig. 15 
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there are several septal stable intervals. In the symmetry requirement we 
prefer choosing the optimal sifting times in the stable intervals, though there 
may be some lower sifting times which accord with lower v. By the way, 
it is not the case that more times of sifting leads to better decomposition. 
On the one hand, as shown in Fig. 9, the additional sifting may lead addi- 
tional error to the decomposition; On the other hand, as shown in Fig. 15, 
the decomposition may be not uniform convergent about the sifting times. 
Certainly, there is also another case that the stable interval do not appear 
for several hundred times of sifting. This is possibly caused by a too small 
value of e. At this time, we suggest replacing the value of e by a bigger one 
and redoing it. In addition, we note that in the stable interval the decompo- 
sition is insensitive to the sifting times and, on the contrary, in the unstable 
interval it differs much. Especially, when the variance ratio v > 100%, the 
corresponding decomposition may be very bad. 

4-4- Effect of Least Extreme-Point Number to the Decomposition 

To stop the decomposition it requires a least number of extreme 
points. This number, denoted by m R , may influence the shape of R. In or- 
der to satisfy the request of boundary linear interpolation we usually choose 
m-R > 4. The default one is 4. If the so-called optimal R for the default 
case has much difference from the data (reflected by a very high u) we can 
increase m R and redo it to get a better one. Fig.9 and 10 are the default 
decompositions. Fig. 11 and 12 are the decompositions with m R = 6. Fig. 15 
and 16 are the decompositions with m R = 8. In Fig. 16 if itlr is increased 
to 20 the AGM curve may become the combination of Mode5 and R which 
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should be avoided, though the variance ratio is much lowered down at this 
time. 
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Figure 17: The odd-even type of 
extreme-point symmetry for Model of 
Fig. 12 (symmetry between the inner 
curves with * and o). Here two outer 
envelop curves are also shown. 



Figure 18: The decomposition result 
for the temperature data given by 
ESMDJII with 11 sifting times, here the 
horizontal axis stands for the time (day). 



5. Performance of ESMDJII 

ESMDJII does the sifting process by using 3 curves Li,L 2 and L 3 
interpolated by the midpoints numerated by 3k + 1, 3k + 2 and 3(k + 1) 
(k = 0, l,---), respectively For this case the mean curve is defined as 
L* = (L\ + L2 + L3) /3. To make a mode symmetric it requires \L*\ < e. This 
kind of symmetry is an extreme-point symmetry in the more extensive mean- 
ing, which only requires the symmetry between L\ + L2 and L3. Particularly, 
when the component is a periodic function it becomes true extreme-point 
symmetry; when the component is a weighted-periodic function it degener- 
ates to the odd-even type of extreme-point symmetry. So ESMDJII can give 
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almost the same decomposition result as ESMDJI for Example 3. By the 
way, for this case ESMDJII only needs 5 sifting times, which is much less 
than 29 for ESMDJI. 

It follows from Fig. 18 that ESMDJII also outputs good AGM curve for 
the temperature data. Its 11 times sifting yields a variance ratio 36.97%, 
which is slightly smaller than 37.11% obtained by ESMDJI with 29 times 
sifting. By making comparison with Fig. 16 we see ESMDJII yields less 
modes than that of ESMDJI. One reason is that, relative to the odd-even 
interpolation, the 3-curve one leads much quicker decreasing to the number 
of extreme points; Another reason is that, ESMDJII adopts a lower type of 
symmetry. 

6. Direct Interpolating Approach for Instantaneous Frequency and 
Amplitude 

Now that the period should be defined relative to a segment of time 
and the frequency needs to be understood point by point, we conciliate this 
conflict by developing a "direct interpolating (DI)" approach for it. Just as 
analyzed in Section 1.4, the instantaneous frequency should be capable of 
reflecting the intermittent case in Fig J, rather than excluding the adjacent- 
equal situation. This request leads some complexity to the processing. Rela- 
tively, the interpolation method for instantaneous amplitude is much simpler. 
The detailed interpolation algorithm is as follows. 
6.1. Interpolation Algorithm 

For the decomposed n IMFs we calculate their instantaneous frequen- 
cies by implementing the following algorithm: 
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Step 1: For each IMF (denoted by (t(k),y(k)) with 1 < k < N) find all 
the interpolating points which satisfy: y(k) > y(k — 1) & y(k) > y(k + 1), 
y(k) > y{k - 1) & y(k) > y(k + 1), y(k) < y(k - 1) & y(k) < y{k + 1) or 
y(k) < y(k — 1) & y{k) < y(k + 1) and numerate them by Ei(ti,yi) with 
1 < i < m. 

Step 2: If there are two adjacent such that y^i = yi or yi = y^+i, 
then define the frequency interpolation coordinates by a« = t(i),fi = 0; 
further if and E i+1 are adjacent extreme points, then define Oj_i = (tj + 
ti- 2 )/2, = 1/ (ti - tj„ 2 ) and a i+2 = (^+3 + U + i)/2, f i+2 = 1/ (^+3 - U +1 ); 
else if and E i+1 (or and are not adjacent extreme points, then 
define a;_i = = l/[(t i+2 - ^-2) - (^+1 - ^)] and a i+2 = t i+2 , fi+ 2 = 

l/[(**+3-*i-i)-(*i+i-*»)]; else > define a; = (£ m +^_i)/2, = — ^_i) . 

Step 3: To add the boundary points with linear interpolating method. For 
the left one, if it is an adjacent equal case, assign the value f± — to ai = t(l); 
if not, assign the value f\ = (f 3 - / 2 )(ai - a 2 )/ (a 3 - a 2 ) + f 2 to a : = t(l); 
further if /1 < then assign /1 = 1/ [2(t 2 — ii)] to a± = t(l). For the right 
one, if it is an adjacent equal case, assign the value f m = to a m = t(N); if 
not, assign the value f m = (/ m -i -/ m - 2 )(a m -a m -i)/(a ) + fm-l tO 

a m = t(N); further if f m < then assign f m = l/[2(i m -i m _i)] to a rn = t(N). 

Step 4: To make the interpolation with all the discrete points (a i; fi) and 
get a curve f(t). To be meaningful, we define the instantaneous frequency 
curve by max{0, fit)}. 

Step 5: To output subplot frequency figures for all IMFs. 
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In addition, since the interpolating method for instantaneous amplitude 
is much simpler we omit its algorithm here. For an IMF, its instantaneous 
amplitude curve can be understood as the upper envelop interpolated by all 
the maxima points of this IMF in the absolute-value form (rely on all the 
extreme points of the original one). In fact, for an IMF obtained under the 
envelop-symmetry scheme or the odd-even extreme-point symmetry scheme, 
its amplitude curve is almost equivalent to the upper envelop of the IMF 
itself with slow modulation. But for the three-curve type the amplitude may 
have too quick modulation which is not preferred. In addition, to reflect the 
energy variation in a distinct way, the instantaneous amplitude and frequency 
can be figured out together. 

6. 2. Performance of DI Approach 

In the following we test the DI approach with the decomposition result 
given in Fig. 12. By implementing the previous algorithm on the 4 IMFs it 
yields an instantaneous frequency and amplitude distribution in Fig. 19. It 
follows from Fl that there are several segments on which the instantaneous 
frequency attains the Nyquist frequency fx = //2, here / = 20Hz stands for 
the sampling rate of the wind data. For a given IMF, such as the third one, 
this kind of figure can reflect clear variation of the frequency and amplitude. 
The comparison between F3 and A3 shows that, at t — 10s, the third IMF 
has a sharp oscillation with a very small amplitude. In addition, the Fig. 20 
shows an intuitional frequency distribution for the 4 IMFs. Notice that the 
frequency-overlapping phenomenon only occurs on the first and second IMFs 
at the time 2.5s, we can accept them as different modes. 
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Figure 19: The instantaneous frequency Figure 20: The frequency distribution of 

(F) and amplitude (A) variations of the the IMFs for the wind data. 

IMFs for the wind data along the time. 
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Figure 21: The time- variation of the to- 
tal energy for all the IMFs of the wind 
data. 
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Figure 22: The variation of the optimal 
AGM curve for the wind data along the 
time. 
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Since the frequency and amplitude (energy) of each IMF shift at any time, 
it is unreasonable for projecting the total energy onto each fixed frequency 
as a Fourier frequency-spectrum or a Hilbert time-frequency-spectrum, after 
all, the total energy of all the IMFs itself changes along the time. With 
this understanding we abandon the spectrum method and turn to discussing 
the time-variation of the total energy. According to the state in Section 
1.3, the j'-th IMF almost accord with the mathematical expression Xj(t) = 
Aj(t) cos 6 j(t) (1 < j < n), here Aj(t) is actually the called amplitude curve. 
Take ESMDJI as a default decomposition, then the odd-even extreme-point 
symmetry scheme assigns Aj(t) a slow modulation feather. Based on this 
consideration, we define the saying total energy in kinetic energy form: 



Certainly, here the word "total energy" is in a generalized sense. For a 
temperature data, it can be understood as the whole oscillation intensity of 
the temperature. By using this definition, we get the corresponding variation 
of the total energy for the wind data. It follows from Fig. 21 that the total 
energy of IMFs has three large peaks in the 15s time segment. To make 
comparison with Fig.22 we see these peaks accord just right with the sunk 
parts of the adaptive global mean. This is a very interesting phenomenon. 
Perhaps it is caused by the energy transfer; perhaps it is just a coincidence; 
perhaps there are other deep causes. This topic is also an attractive one to 
be discussed. 




(8) 



i=i 
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7. Results 

Here we have proposed an "extreme-point symmetric mode decompo- 
sition (ESMD)" method for nonlinear and non-stationary signal processing. 
There are two parts for it. The first part is the decomposition approach 
which yields a series of intrinsic mode functions (IMFs) together with an 
optimal "adaptive global mean (AGM)" curve, the second part is the "di- 
rect interpolating (DI)" approach which yields instantaneous amplitudes and 
frequencies for the IMFs together with a time- varying figure for the total en- 
ergy. It can be seen as a new alternate of the well-known "Hilbert-Huang 
transform (HHT)" method. Relative to the HHT method it has five charac- 
teristics as follows: 

(1) Different from constructing 2 outer envelopes, its sifting process is imple- 
mented by the aid of 1, 2, 3 or more inner curves interpolated by the mid- 
points of the line segments connecting the local maxima and minima points. 
Accordingly, the ESMD is classified into ESMDJ, ESMD JI, ESMDJII, etc. 

(2) It does not decompose the signal to the last trend curve with at most one 
extreme point, it optimizes the residual component to be an optimal AGM 
curve which possesses a certain number of extreme points. By means of this 
optimizing process one can determine the sifting times and output the opti- 
mal decomposition. 

(3) The concept of symmetry for ESMD has wider extension than the en- 
velop symmetry. From the viewpoint of material movement, as a matter of 
fact, the oscillation occurs around the equilibrium which may also shift its 
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location during this process. So the extreme-point symmetry actually re- 
flects the local symmetry about itself. ESMDJ adopts a rigid extreme-point 
symmetry which is more rigorous than envelop symmetry; ESMD_II adopts 
an odd-even type of extreme-point symmetry which is almost equivalent to 
the envelop symmetry; ESMDJII adopts a three-curve type of extreme-point 
symmetry which is more general than the envelop symmetry. 

(4) The definition of IMF is extended. The new form not only includes the 
intermittent case but also liberalizes the request on symmetry. 

(5) The Hilbert-spectral-analysis approach for instantaneous frequency and 
amplitude is substituted by the data-based DI one. This new approach easily 
conciliates the conflict: the period should be defined relative to a segment of 
time and the frequency needs to be understood point by point. It is unreason- 
able for projecting the total energy onto each fixed frequency as a Fourier 
frequency-spectrum or a Hilbert time-frequency-spectrum, after all, the total 
energy of all the IMFs itself changes along the time. The DI approach can 
not only yield clear distribution for instantaneous frequency and amplitude 
but also reflect the time-variation of total energy in a distinct way. 

By making comparison among ESMDJ, ESMDJI and ESMDJII we get 
some understandings on the effect of interpolation, that is, as the number of 
interpolating curves increases, 1) the mode number decreases along; 2) the 
symmetry degree decreases along; 3) the amplitude modulation increases 
along; 4) the decomposition efficiency increases along (needs less sifting 
times). With these understandings, we prefer doing the decomposition with 
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the eclectic one. In fact, ESMDJ can only yield acceptable decomposition 
with imperfect sifting at low times. The equal-amplitude phenomenon will 
occur at high sifting times which should be avoided in consideration of phys- 
ical meaning. Though ESMDJII has high decomposition efficiency, its low 
degree of symmetry and quick amplitude modulation are disadvantageous for 
frequency and energy analysis. The signal processing tests also indicate that 
ESMDJI is superior to ESMDJ and ESMDJII and its decomposition result 
is more preferable. 

Though these three types of interpolation have much difference, they have 
a mutual advantage. Not only ESMDJI but also ESMDJ and ESMDJII can 
output good AGM curve. In fact, this advantage owes to the feature of inner 
interpolation. Besides the decomposition, the ESMD method also offers a 
good adaptive approach for finding the optimal global mean fitting curve to a 
given data. It is superior to the commonly used least-square method and the 
running mean approach. In fact, the first one is awkward in application due 
to the request on a priori function form, the second one lacks of theoretical 
base and different choices of weight coefficients may result in different curves 
with two uncertain boundaries. 

We note here that the type of interpolation is not an essential factor, the 
approach developed here is also suitable for the envelop-symmetric scheme 
adopted by the EMD method [see Appendix B]. From the viewpoint of stop- 
page criterion our strategy is an synthetic one. We have chosen a preset- 
error condition to ensure the symmetry and an ensemble optimal-sifting- 
times (OST) approach to optimize the whole decomposition. In addition 
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to this stoppage criterion there is also another applicable one. In fact, in 
view of the intermittent feature of IMF's symmetry [Wang and Li (2012)], 
we can abandon the preset-error condition and implement the OST approach 
repeatedly and draw out a series of IMFs. If the mode symmetry is the sole 
attention, this approach is a better choice. But it follows from the test in 
Appendix C that the last residual may be an inferior one. To make up this 
defect, one can further implement the ensemble OST approach on the whole 
IMF sets. Certainly, its time-consuming would be longer than the present 
one. 

Acknowledgments. We would like to thank Professor Norden E. Huang 
for his enthusiastic encouragement and many stimulating discussions on the 
topics related to the present research. 

Appendix A: Boundary Processing 

In our program codes we have developed the linear interpolation method 
given by Wu and Huang (2009) and revised the interpolation styles for the 
too steep boundary case. This revision can make the boundary much more 
stable, even for the tests given by Wang and Li (2012) with 100, 000 sifting 
times. Take the left boundary processing as an example. Let y(t) = kit + b\ 
and y(t) = k 2 t + b 2 be the upper and lower lines interpolated by the first two 
maxima and minima points, respectively. Also denote the first point of the 
data by Y\. According to Wu and Huang's classification, (1) if b 2 < Yi < 
bi, then define b\ and b 2 as the boundary maximum and minimum points, 
respectively; (2) if Y\ > b\ (or Y\ < b 2 ), then define Y\ and b 2 (or b\ and Y\) 
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as the boundary maximum and minimum points. But if the boundary is too 
steep (relative to these two interpolating lines), this kind of processing may 
lead instability to the decomposition. So we substitute the second term by: 
(2) if h <Y 1 < (61 + 62) /2 + (61 - 6 2 ) = (36 x - 6 2 )/2 (or (3b 2 - b x )/2 = 
(61 + b 2 )/2 — (61 — b 2 ) < Y\ < b 2 ), then define Yi and b 2 (or bi and Yi) as the 
boundary maximum and minimum points; (3) if Y\ > (3bi — b 2 )/2 (or Y\ < 
(3b 2 — bi)/2), then define Y\ as the boundary maximum (minimum) point and 
take the boundary minimum (maximum) point by using new interpolating 
line from the first minimum (maximum) point: y(t) = k*t + b*, here k* relies 
on the point (0, Y{) and the first maximum (minimum) point. The detailed 
processing are depicted in the following figure. 
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Figure 23: The developed linear interpolation method for the left boundary processing. 
Case A: b 2 <Y 1 < bi\ Case B: b 1 <Y 1 < (36 x - b 2 )/2; Case C: Y 1 > (36 a - b 2 )/2. 

Appendix B: Decomposition with Envelop-Symmetric Scheme 
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The approach developed in this paper is also suitable for the envelop- 
symmetric scheme adopted by the EMD method. At this time, the decom- 
position result is similar to that of odd-even extreme-point symmetric one. 
To make comparison with ESMDJI at 30 sifting times, we choose the second 
optimal one 39 in the interval [1, 100] (the first optimal one is 62). For this 
case the AGM curve is as good as ESMDJI with variance ratio v = 33.8% 
and Model and Mode2 are very similar to that of Fig. 12. The difference lies 
in the low- frequency modes (Mode3 and Mode4). This indicates the whole 
outer symmetry is similar to the local inner symmetry for the case with dense 
extreme points. But for the sparse case, the outer and inner interpolations 
may result in different results. Notice that the magnitudes of the upper and 
lower outer envelopes are bigger than that of odd and even inner curves, their 
interpolating uncertainty should be bigger than the inner ones for the sparse 
case. Hence, the result given by ESMDJI should be more credible. 

Appendix C: Decomposition with Optimal-Sifting-Times Approach 

In the following we try another stoppage criterion which is associated 
with the single usage of optimal-sifting-times (OST) approach. For a preset 
maximum sifting times K max , we can select an optimal one in the integer 
interval [l,A" maa; ] which accords with the minimum value of A max , where 
A m ax stands for the maximum amplitude of an quasi- mode's mean cure L*. 
By implementing this OST approach repeatedly we can draw out a series of 
IMFs. For convenience of processing we take the envelop-symmetric scheme 
as an example. At this time, L* is simply the mean of the upper and lower 
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Figure 24: Under the envelop-symmetric 
scheme, the distribution of the variance 
ratio v = a/ao along the sifting times for 
the wind data. 



Figure 25: Under the envelop-symmetric 
scheme, the decomposition result for the 
wind data with 39 sifting times, here the 
horizontal axis stands for the time (sec- 
ond). 



envelops and the corresponding decomposition results are given in Fig. 26 and 
27 with K max = 100 and 200, respectively. Though in the interval [1, 200] we 
have more choice and the symmetry of IMFs should be better, its last residual 
R is worse than that in [1, 100]. That is to say, the optimal processing on 
IMFs can not ensure an optimal global mean cure. This defect can be made 
up by further optimizing K max with respect to the variance ratio v = cr/exo- 
Certainly, it would cost a longer time-consuming than the one adopted by 
this paper due to more times of calculation. 

We also note that the degree of symmetry does not increase uniformly 
along the sifting times. Just as indicated by Wang and Li (2012), its variation 
behaves in an intermittent manner [see Fig. 28]. 
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Figure 26: Under the envelop-symmetric 
scheme, the decomposition result for the 
wind data given by optimal-sifting-times 
(OST) approach with K max = 100, here 
the horizontal axis stands for the time 
(second) . 



Figure 27: Under the envelop-symmetric 
scheme, the decomposition result for the 
wind data given by optimal-sifting-times 
(OST) approach with K max — 200, here 
the horizontal axis stands for the time 
(second). 




Sifting Times 



Figure 28: The variation of maximum amplitude for the mean curve of IMF1 along the 
sifting times, here the wind data is used. 
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